#include "fuel_cell_infos.h"

#include <stdio.h>
#include <math.h>
#include "point.h"
//~ #include "mpi_tools.h"
#include "opencl_tools.h"
#include "matrix.h"
#include "physical_constants.h"
#include "params.h"
#include "list.h"

Point3 nbNodes;
Point3 matrixSize;

int   ***fuelCellMap = NULL;
tData ***heatCapacity = NULL;

ListInt *packedListX = NULL;
ListInt *packedListY = NULL;
ListInt *packedListZ = NULL;
int *ePackedListX = NULL;
int *ePackedListY = NULL;
int *ePackedListZ = NULL;
ListInt *reversePackedListX = NULL;
ListInt *reversePackedListY = NULL;
ListInt *reversePackedListZ = NULL;
Point3 packedListSize;

CoeffBlock ***blocksCoeff = NULL;

void InitFuelCell();

void InitBlocksCoeff();
void ComputeBlockCoeff(Point3 packedCoord);
Point3 PackedCoordToBlockCoord(Point3 packedCoord);
// Transfert
// void ComputeGazInfos(Point3 packedCoord, Point3 blockCoord, CoeffBlock *cb);
tData ComputeTransfert(int mat1, int mat2, int mat3, tData surface, tData d1, tData d2, tData d3, Point3 blockCoord, tData heatCapacity, int direction);
tData ComputeTransfertBlock(Point3 blockCoord1, Point3 blockCoord2, Point3 blockCoord3, int direction);
tData ComputeConstant(Point3 blockCoord);
tData ComputeDebit(int mat);
// Material
int GetMaterialOf(Point3 blockCoord);
// Dimensions
void ComputeGlobalMatrixSize();
tData GetSizeX(int blockPosX);
tData GetSizeY(int blockPosY);
tData GetSizeZ(int blockPosZ);
// Heat
void InitHeatCapacity();
tData GetHeatCapacity(Point3 blockCoord);
tData ComputeHeatCapacity(int x, int y, int z);
tData getVgaz(int mat);
tData getDebit_times_T_gaz(int mat, tData surface);
// tData GetM_div_T_gaz(int mat);

void InitFuelCellInfos(int *granularity){
    InitFuelCell(granularity);
    InitHeatCapacity();
    InitBlocksCoeff();
}

void FreeBlocksCoeff(){
    if( packedListX != NULL )
        free(packedListX);
    if( packedListY != NULL )
        free(packedListY);
    if( packedListZ != NULL )
        free(packedListZ);

    if( reversePackedListX != NULL )
        free(reversePackedListX);
    if( reversePackedListY != NULL )
        free(reversePackedListY);
    if( reversePackedListZ != NULL )
        free(reversePackedListZ);

    if( ePackedListX != NULL )
        free(ePackedListX);
    if( ePackedListY != NULL )
        free(ePackedListY);
    if( ePackedListZ != NULL )
        free(ePackedListZ);

    packedListX = NULL;
    packedListY = NULL;
    packedListZ = NULL;
    ePackedListX = NULL;
    ePackedListY = NULL;
    ePackedListZ = NULL;
    reversePackedListX = NULL;
    reversePackedListY = NULL;
    reversePackedListZ = NULL;
}

void InitBlocksCoeff(){
    FreeBlocksCoeff();

    packedListX = ListBuildPacked( blockNbNodesX );
    packedListY = ListBuildPacked( blockNbNodesY );
    packedListZ = ListBuildPacked( blockNbNodesZ );

    reversePackedListX = ListBuildReversePacked( blockNbNodesX );
    reversePackedListY = ListBuildReversePacked( blockNbNodesY );
    reversePackedListZ = ListBuildReversePacked( blockNbNodesZ );

    ePackedListX = ListExpand( packedListX );
    ePackedListY = ListExpand( packedListY );
    ePackedListZ = ListExpand( packedListZ );

    packedListSize = Point3i(
            packedListX->size - 1, packedListY->size - 1, packedListZ->size - 1);

    blocksCoeff = (CoeffBlock ***)Matrix3Allocate(
            Point3i(1, 1, 1), packedListSize, sizeof(CoeffBlock) );

    for(int x = 1; x < packedListSize.x; x++)
        for(int y = 1; y < packedListSize.y; y++)
            for(int z = 1; z < packedListSize.z; z++)
                ComputeBlockCoeff( Point3i(x, y, z) );
}

void ComputeBlockCoeff(Point3 packedCoord){
    Point3 blockCoord = PackedCoordToBlockCoord(packedCoord);
    tData heatCapacity = GetHeatCapacity(blockCoord);
    CoeffBlock *cb = &blocksCoeff[packedCoord.x][packedCoord.y][packedCoord.z];
    Point3 blockCoord2, blockCoord3;

    int mat = GetMaterialOf(blockCoord);
    int matType = tMaterial[mat].type;
    cb->constant = 0;
    cb->special_type = 0;
    if(mat == MAT_FICTIF){
      cb->self = 0;
      tData coef = -1;
      for(int i=0; i<6; i++){
            blockCoord2 = PackedCoordToBlockCoord( Point3Add(packedCoord, direction3[i]) );
            blockCoord3 = PackedCoordToBlockCoord( Point3Sub(packedCoord, direction3[i]) );
            cb->transfert[i] = ComputeTransfertBlock( blockCoord, blockCoord2, blockCoord3, i );

            int mat2 = GetMaterialOf(blockCoord2);
            int mat3 = GetMaterialOf(blockCoord3);
            if((mat2 == mat3 && mat2 != MAT_FICTIF) || (mat2 == MAT_COOLER && mat3 != MAT_OUTSIDE) || (mat3 == MAT_COOLER && mat2 != MAT_OUTSIDE)) {
              if(mat2 == MAT_H2 || mat2 == MAT_AIR){
                cb->special_type = 1;
              }
              coef = cb->transfert[i];
              if(coef == 0)
                die_cl("Error, coef can't be 0 here, mat2: %d, mat3: %d", mat2, mat3);
              // This is not an interface, it should have a constant like both its neighbours
              cb->constant = deltaT * ComputeConstant(blockCoord2) / heatCapacity;
            }
       }
       if(coef != -1){
          // this point is fictif but is not an interface. So we fix the transferts
         for(int i=0; i<6; ++i){
           if(cb->transfert[i] == 0)
            cb->transfert[i] = coef;
         }
         if(cb->special_type == 0){
           cb->self = 1;
           for(int i=0; i<6; i++){
              cb->self -= cb->transfert[i];
           }
         }
       }
//        printf("cbself: %e\n", cb->self);
    }else{
      if(matType == TYPE_GAZ){
        cb->special_type = 1;
        cb->constant = 0;
        if(mat == MAT_OUTSIDE)
          die_cl("Error: Mat outside spotted\n");
        for(int i=0; i<6; i++){
          if(i%2 != 0){
            blockCoord2 = PackedCoordToBlockCoord( Point3Add(packedCoord, direction3[i]) );
            blockCoord3 = PackedCoordToBlockCoord( Point3Sub(packedCoord, direction3[i]) );
            cb->transfert[i] = ComputeTransfertBlock( blockCoord, blockCoord2, blockCoord3, i );
//             cb->self -= cb->transfert[i];
          }
        }
      }else{
        // Liquid or Solid
        cb->self = 1;
        for(int i=0; i<6; i++){
            blockCoord2 = PackedCoordToBlockCoord( Point3Add(packedCoord, direction3[i]) );
            blockCoord3 = PackedCoordToBlockCoord( Point3Sub(packedCoord, direction3[i]) );
            cb->transfert[i] = ComputeTransfertBlock( blockCoord, blockCoord2, blockCoord3, i );
            if(! (cb->transfert[i] >= 0 && cb->transfert[i] < 1000))
              die_cl("invalid transfert value: %e, mat1: %d, mat2: %d, mat3: %d\n", cb->transfert[i], mat, GetMaterialOf(blockCoord2), GetMaterialOf(blockCoord3));
            cb->self -= cb->transfert[i];
        }
        // Constant (only for reactAnode & reactCathode)
        cb->constant = deltaT * ComputeConstant( blockCoord ) / heatCapacity;
//         if(cb->constant != 0)
//           printf("constante: %e\n", cb->constant/deltaT);
      }
    }
}

tData getVgaz(int mat){
  tData Pr, Re;
  if(mat == MAT_AIR){
    Pr = PHYS_nuV_AIR * tMaterial[mat].CP / tMaterial[mat].K;
    Re = pow(PHYS_Nu_gaz / (0.023*pow(Pr, 0.4)), 1./0.8);
//     printf("Re_Air= %e\n", Re);
//     printf("V_Air = %e\n", Re * PHYS_nuV_AIR / PHYS_Dh);
    return 1.2;
    return Re * PHYS_nuV_AIR / PHYS_Dh;
  }
  if(mat == MAT_H2){
    Pr = PHYS_nuV_H2 * tMaterial[mat].CP / tMaterial[mat].K;
    Re = pow(PHYS_Nu_gaz / (0.023*pow(Pr, 0.4)), 1./0.8);
//     printf("Re_H2= %e\n", Re);
//     printf("V_H2 = %e\n", Re * PHYS_nuV_H2 / PHYS_Dh);
    return 1.2;
    return Re * PHYS_nuV_H2 / PHYS_Dh;
  }
  die_cl("Unknown gaz: %d", mat);
  return 999;
}


tData getDebit_times_T_gaz(int mat, tData surface){
    if( mat == MAT_H2 )
      return tMaterial[mat].Rho * surface * getVgaz(mat);
//         return  (PHYS_P_A * PHYS_M_H2 * getVgaz(mat)) / PHYS_R;
    if( mat == MAT_AIR )
      return tMaterial[mat].Rho * surface * getVgaz(mat);
//         return  (PHYS_P_C * PHYS_M_AIR * getVgaz(mat)) / PHYS_R;

    die_cl("getDebit_times_T_gaz mat = %d", mat);
    return 666;
}

//     tData GetM_div_T_gaz(int mat){
//         if( mat == MAT_H2 )
//             return PHYS_R / PHYS_P_A / PHYS_UT_H2 / ( 1 + PHYS_HUMIDITY_H2 ) * PHYS_I / (2 * PHYS_F) / 2 / NBNODES_Z;
//         if( mat == MAT_AIR )
//             return PHYS_R / PHYS_P_C / (PHYS_UT_AIR * PHYS_X_O2) / ( 1 + PHYS_HUMIDITY_AIR ) * PHYS_I / (4 * PHYS_F) / 2 / NBNODES_Z;
//
//         die_cl("GetM_div_T_gaz mat = %d", mat);
//         return 666;
//     }

// void ComputeGazInfos(Point3 packedCoord, Point3 blockCoord, CoeffBlock *cb){
//     int mat = GetMaterialOf(blockCoord);
//     tData heatCapacity = GetHeatCapacity(blockCoord);
//     tData cp = tMaterial[mat].CP;
//     cb->self = 0;
//     cb->special_self= 0;
//
//     for(int i=0; i<6; i++){
//         Point3 blockCoord2 = PackedCoordToBlockCoord( Point3Add(packedCoord, direction3[i]) );
//         Point3 blockCoord3 = PackedCoordToBlockCoord( Point3Sub(packedCoord, direction3[i]) );
//         int mat2 = GetMaterialOf(blockCoord2);
//
//         if( mat == mat2 ){
//             cb->transfert[i] = 0.;
//             cb->special[i] = gazOrientation[i] * deltaT * cp * GetM_div_T_gaz(mat) / heatCapacity;
//             cb->special_self -= cb->special[i];
//         }else{
//             cb->transfert[i] = deltaT * ComputeTransfertBlock(blockCoord, blockCoord2, blockCoord3, i) /
//                 heatCapacity;
//             cb->special[i] = 0.;
//             cb->self -= cb->transfert[i];
//         }
//     }
// }

tData ComputeTransfertBlock(Point3 blockCoord1, Point3 blockCoord2, Point3 blockCoord3, int direction){
    int mat1 = GetMaterialOf(blockCoord1);
    int mat2 = GetMaterialOf(blockCoord2);
    int mat3 = GetMaterialOf(blockCoord3);
    tData heatCapacity = GetHeatCapacity(blockCoord1);
    Point3 dir = Point3Abs(direction3[direction]);
    tData surface = 0;

    tData d1 =  \
                GetSizeX(blockCoord1.x) * dir.x +
                GetSizeY(blockCoord1.y) * dir.y +
                GetSizeZ(blockCoord1.z) * dir.z;
    tData d2 =  \
                GetSizeX(blockCoord2.x) * dir.x +
                GetSizeY(blockCoord2.y) * dir.y +
                GetSizeZ(blockCoord2.z) * dir.z;

    tData d3 =  \
                GetSizeX(blockCoord3.x) * dir.x +
                GetSizeY(blockCoord3.y) * dir.y +
                GetSizeZ(blockCoord3.z) * dir.z;

      if( dir.x == 1 )
        surface =
            GetSizeY(blockCoord1.y) *
            GetSizeZ(blockCoord1.z);
      else if( dir.y == 1 )
          surface =
              GetSizeX(blockCoord1.x) *
              GetSizeZ(blockCoord1.z);
      else if( dir.z == 1 )
          surface =
              GetSizeX(blockCoord1.x) *
              GetSizeY(blockCoord1.y);
      else
          die_cl(" Error surface.");
    if(surface == 0){
      if( dir.x == 1 )
          surface =
              GetSizeY(blockCoord2.y) *
              GetSizeZ(blockCoord2.z);
      else if( dir.y == 1 )
          surface =
              GetSizeX(blockCoord2.x) *
              GetSizeZ(blockCoord2.z);
      else if( dir.z == 1 )
          surface =
              GetSizeX(blockCoord2.x) *
              GetSizeY(blockCoord2.y);
    }else{
//       printf("x: %e, y: %e\n", GetSizeX(blockCoord1.x), GetSizeX(blockCoord1.y));
      return ComputeTransfert(mat1, mat2, mat3, surface, d1, d2, d3, blockCoord1, heatCapacity, direction);
    }
    if(surface == 0){
      if( dir.x == 1 )
          surface =
              GetSizeY(blockCoord3.y) *
              GetSizeZ(blockCoord3.z);
      else if( dir.y == 1 )
          surface =
              GetSizeX(blockCoord3.x) *
              GetSizeZ(blockCoord3.z);
      else if( dir.z == 1 )
          surface =
              GetSizeX(blockCoord3.x) *
              GetSizeY(blockCoord3.y);
    }else{
      return ComputeTransfert(mat1, mat2, mat3, surface, d1, d2, d3, blockCoord2, heatCapacity, direction);
    }
    if(surface == 0)
      surface = 2.953271*pow(10, -5);
    return ComputeTransfert(mat1, mat2, mat3, surface, d1, d2, d3, blockCoord3, heatCapacity, direction);
}

tData ComputeTransfert(int mat1, int mat2, int mat3, tData surface, tData d1, tData d2, tData d3, Point3 blockCoord, tData heatCapacity, int direction){
    int matType1 = tMaterial[mat1].type;
    int matType2 = tMaterial[mat2].type;
    int matType3 = tMaterial[mat3].type;
    tData deltaD = (d1+d2)/2.;
    if(deltaD == 0)
      deltaD = d3/2.;
    tData deltaD2 = deltaD * deltaD;
//     Point3 dir = Point3Abs(direction3[direction]);
//     printf("Computing transfert\n");
    if(blockCoord.x == 0 && blockCoord.y == 2 && blockCoord.z == 5){
      printf("mat1:%d, mat2: %d, mat3: %d\n", mat1, mat2, mat3);
    }

    if( mat1 == MAT_OUTSIDE)
      die_cl("Error, souldn't have outside here");

    // Pas de transfert entre matériaux fictifs.
    if( mat1 == MAT_FICTIF){
      if(mat2 == MAT_FICTIF || mat3 == MAT_FICTIF ){
        return 0;
      }
    }

    //
    // Interface
    //
    if(mat1 == MAT_FICTIF && mat2 != mat3){
      if(mat2 == MAT_OUTSIDE){
        // Solide -> exterieur
        tData K = tMaterial[mat3].K;
        if(K == 0)
          die_cl("K ne peut pas valoir 0 ici, mat3 vaut: %d", mat3);
//         printf("solide -> ext: %e\n", (deltaD * tMaterial[mat2].H * S) / (K + tMaterial[mat2].H * S * deltaD));
        return (deltaD * tMaterial[mat2].H * surface) / (K + tMaterial[mat2].H * surface * deltaD);
      }
      if(matType2 == TYPE_GAZ){
        // gaz -> solide
        tData K = tMaterial[mat3].K;
        //printf("gaz -> solide: %f\n", (tMaterial[mat2].H * S) / ((tMaterial[mat2].H * S) + (K / deltaD)));
        return (tMaterial[mat2].H * surface) / ((tMaterial[mat2].H * surface) + (K / deltaD));
      }
      // mat2 est solide
      if(mat3 == MAT_OUTSIDE){
        // exterieur -> solide
        tData K = tMaterial[mat2].K;
        if(K == 0)
          die_cl("K ne peut pas valoir 0 ici, mat2 vaut: %d", mat2);
//         printf("ext -> solide: %f\n", K / (K + deltaD * tMaterial[mat3].H * S));
        return K / (K + deltaD * tMaterial[mat3].H * surface);
      }
      // Interface Solide - liquide
      if(matType2 == TYPE_LIQUID){
//         printf("Solide - liquid: %e, surface: %e\n", surface * tMaterial[mat2].H * ALPHA_LIQUID, surface);
        return surface * tMaterial[mat2].H * ALPHA_LIQUID;
      }
      if(matType3 == TYPE_LIQUID){
//         printf("Solide - liquid: %e, surface: %e\n", surface * tMaterial[mat3].H * ALPHA_LIQUID, surface);
        return surface * tMaterial[mat3].H * ALPHA_LIQUID;
      }
      if(matType3 == TYPE_GAZ){
        // solide -> gaz
        tData K = tMaterial[mat2].K;
        //printf("solide -> gaz: %f", K / ((tMaterial[mat3].H * S * deltaD) + K));
        return K / ((tMaterial[mat3].H * surface * deltaD) + K);
      }
      // Interface Solide1 <- fictif ->  Solide 2
      tData K2 = tMaterial[mat2].K;
      tData K3 = tMaterial[mat3].K;
      //printf("solide1 - solide2: %f\n", (K2*d3)/(d3*K2+d2*K3));
      return (K2*d3)/(d3*K2+d2*K3);
    }
    //
    // Pas une interface
    //
    if( mat1 == MAT_FICTIF){
      // Bloc fictif entre deux blocs identiques
      // Traiter comme les blocs voisins
      // Info: mat2 == mat3
      mat1 = mat2;
      matType1 = matType2;
    }

    if( matType1 == TYPE_LIQUID){
      // Liquide -> Liquide
//       printf("Liquid - liquid: %e\n", ALPHA_LIQUID * tMaterial[mat1].CP * (GetSizeX(blockCoord.x) + GetSizeY(blockCoord.y)) * PHYS_M_EAU / NBNODES_Z);
      if(GetSizeX(blockCoord.x) == 0 && GetSizeY(blockCoord.y) == 0){
//         printf("(sp) liquid - liquid: %e\n", ALPHA_LIQUID * tMaterial[mat1].CP * 1.38 * pow(10, -5) * PHYS_M_EAU / NBNODES_Z);
        return ALPHA_LIQUID * tMaterial[mat1].CP * 1.38 * pow(10, -5) * PHYS_M_EAU / NBNODES_Z;
      }
//       printf("liquid - liquid: %e, x: %e, y: %e\n", ALPHA_LIQUID * tMaterial[mat1].CP * (GetSizeX(blockCoord.x) + GetSizeY(blockCoord.y)) * PHYS_M_EAU / NBNODES_Z, GetSizeX(blockCoord.x), GetSizeY(blockCoord.y));
      return ALPHA_LIQUID * tMaterial[mat1].CP * (GetSizeX(blockCoord.x) + GetSizeY(blockCoord.y)) * PHYS_M_EAU / NBNODES_Z;
    }
    if( matType1 == TYPE_GAZ ){
      // Gaz -> Gaz
      //printf("gaz ->gaz: %f\n", getDebit_times_T_gaz(mat1) * tMaterial[mat1].CP * deltaT / deltaD);
      return getDebit_times_T_gaz(mat1, surface) * tMaterial[mat1].CP * deltaT / deltaD;
    }
    // Solide -> Solide
    tData K = tMaterial[mat1].K;
//      if(mat1 == MAT_CATHODE_REACT)
//       printf("solide -> solide: %e\n", (deltaT * K) / (deltaD2 * heatCapacity));
    return (deltaT * K) / (deltaD2 * heatCapacity);
}

tData ComputeConstant(Point3 blockCoord){
    int mat = GetMaterialOf(blockCoord);
//     tData X = GetSizeX(blockCoord.x);
    tData Y = GetSizeY(blockCoord.y);
    tData Z = GetSizeZ(blockCoord.z);

    tData result = 0;

    if( mat == MAT_ANODE_REACT ){
        result = (PHYS_LV * PHYS_ALPHA * PHYS_I/* * Y * Z*/) / (2 * PHYS_F /** X */* Y * Z);
//         printf("react anode : %f", (PHYS_LV * PHYS_ALPHA * PHYS_I/* * Y * Z*/) / (2 * PHYS_F));
    }

    if( mat == MAT_CATHODE_REACT ){
        result = /*Y * Z * */(
                ((PHYS_LV * (2 * PHYS_ALPHA + 1) - PHYS_TDELTAS) / ( 2 * PHYS_F ) * PHYS_I
                + PHYS_B_O2 * PHYS_LOG)) / (/*X**/Y*Z);
//         printf("react cathode : %f", (PHYS_LV * (2 * PHYS_ALPHA + 1) - PHYS_TDELTAS) / ( 2 * PHYS_F ) * PHYS_I
//                 + PHYS_B_O2 * PHYS_LOG);
    }

    return result / blockNbNodesX->list[blockCoord.x];
}

Point3 PackedCoordToBlockCoord(Point3 packedCoord){
    return Point3i(
            reversePackedListX->list[packedCoord.x],
            reversePackedListY->list[packedCoord.y],
            reversePackedListZ->list[packedCoord.z]);
}

void InitHeatCapacity(){
    heatCapacity = (tData ***)Matrix3Allocate(zero3, nbNodes, sizeof(tData));

    for(int x = 0; x < nbNodes.x; x++)
        for(int y = 0; y < nbNodes.y; y++)
            for(int z = 0; z < nbNodes.z; z++)
                heatCapacity[x][y][z] = ComputeHeatCapacity(x, y, z);
}

tData ComputeHeatCapacity(int x, int y, int z){
    int mat = fuelCellMap[x][y][z];
    tData rho = tMaterial[mat].Rho;
    tData cp = tMaterial[mat].CP;
    tData X = GetSizeX(x);
    tData Y = GetSizeY(y);
    tData Z = GetSizeZ(z);
    tData hc;

    hc = cp * rho/* * X * Y * Z*/;

    if( hc == 0.0 )
        die_cl("HeatCapacity too small. coord=(%d, %d, %d). "
                " size( %f, %f, %f) , "
                "cp = %f, rho = %f, Material = %s",
                x, y, z, X, Y, Z, cp, rho, tMaterial[mat].name);

    return hc;
}

tData GetHeatCapacity(Point3 blockCoord){
    return heatCapacity[blockCoord.x][blockCoord.y][blockCoord.z];
}

int GetMaterialOf(Point3 blockCoord){
    return fuelCellMap[blockCoord.x][blockCoord.y][blockCoord.z];
}

tData GetSizeX(int blockPosX){
    return blockSizeX[blockPosX] / blockNbNodesX->list[blockPosX];
}

tData GetSizeY(int blockPosY){
    return blockSizeY[blockPosY] / blockNbNodesY->list[blockPosY];
}

tData GetSizeZ(int blockPosZ){
    return blockSizeZ[blockPosZ] / blockNbNodesZ->list[blockPosZ];
}

